import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

def c_func(w,m): return np.sign(np.cos(w)) * np.abs(np.cos(w))**m
def s_func(w,m): return np.sign(np.sin(w)) * np.abs(np.sin(w))**m

def sample_superellipsoid_parametric(A=1, B=1, C=1, e1=0.5, e2=0.3, u_num_samples=50, v_num_samples=50):
    v = np.linspace(-np.pi/2, np.pi/2, v_num_samples)
    u = np.linspace(-np.pi, np.pi, u_num_samples)
    uu, vv = np.meshgrid(u, v)
    x = A * c_func(vv, e1) * c_func(uu, e2)
    y = B * c_func(vv, e1) * s_func(uu, e2)
    z = C * s_func(vv, e1)
    return np.array([x, y, z]).reshape((3,-1)).T

a,b,c,e1,e2 = 1, 1, 1, 0.5, 0.7
points = sample_superellipsoid_parametric(a,b,c,e1,e2)

fig = plt.figure(figsize=(6,6))
ax = fig.add_subplot(projection='3d')    
ax.set_title("$a={}$, $b={}$, $c={}$, $e_1={}$, $e_2={}$".format(a,b,c,e1,e2))
_=ax.scatter(points[:,0], points[:,1], points[:,2], c=-points[:,2], cmap="gist_rainbow")

fig = plt.figure(figsize=(10,10))
e_values = [1.0, 0.5, 0.01]
for i,e1 in enumerate(e_values):
    for j,e2 in enumerate(e_values):
        ax = fig.add_subplot(len(e_values), len(e_values), i*len(e_values)+j+1, projection='3d')
        points = sample_superellipsoid_parametric(1,1,1,e1,e2)
        ax.set_title("$e_1={}$, $e_2={}$".format(e1,e2))
        _=ax.scatter(points[:,0], points[:,1], points[:, 2], c=-points[:,2], cmap="gist_rainbow")
fig.tight_layout()

# see: https://stackoverflow.com/a/26127012
# Optimized version of this code is provided at the end of this post!
def sample_fibonacci_sphere_(samples): 
    points = np.zeros((samples, 3))
    phi = np.pi * (3 - np.sqrt(5))  # golden angle in radians
    
    for i in range(samples):
        y = 1 - (i / (samples - 1)) * 2   # y goes from 1 to -1
        radius = np.sqrt(1 - y * y)  # radius at y
        theta = phi * i  # golden angle increment
        x = np.cos(theta) * radius
        z = np.sin(theta) * radius
        points[i] = [x, y, z]

    return points

points = sample_fibonacci_sphere_(samples=2500)
fig = plt.figure(figsize=(6,6))
ax = fig.add_subplot(111, projection='3d')
_=ax.scatter(points[:,0], points[:,1], points[:, 2], c=-points[:,2], cmap="gist_rainbow")

# use Solina's method but assume that a=b=c=1
def compute_superellipsoid_distance(x, y, z, e1, e2):
    f = ((np.abs(x))**(2/e2) + (np.abs(y))**(2/e2) )**(e2/e1) + (np.abs(z))**(2/e1)
    d = np.abs(1 - f**(e1/2))
    return np.where(f>1, d, -d) # use sign to find if the point is inside or outside

x = np.linspace(-1.1, 2, 1000)
y = np.linspace(-1.1, 2, 1000)
xx, yy = np.meshgrid(x, y)
zz = 0.0
e_values = [(0.5, 0.5), (0.1, 0.1), (0.01, 0.01)]

fig, axs = plt.subplots(1, len(e_values), figsize=(12,4))
for i, (e1,e2) in enumerate(e_values):
    d = compute_superellipsoid_distance(xx, yy, zz, e1, e2)
    axs[i].contourf(np.abs(d), vmin=0, vmax=1.37) # vmax is found via print(d.min(), d.max())
    axs[i].contour(d, levels=[0], colors="white")
    axs[i].set_title("$e_1={}$, $e_2={}$".format(e1,e2))
fig.suptitle("White pixels represent Superellipsoid surface. Other colors represent the absolute distance to Superellipsoid surface.")
fig.tight_layout()

# Optimized version of this code is provided at the end of this post! Don't use this function!
def sample_superellipsoid_fibonacci_(a, b, c, e1, e2, num_samples=2500, limit_e1_e2=True):
    if limit_e1_e2:
        e1 = np.clip(e1, 0.01, 1.0)
        e2 = np.clip(e2, 0.01, 1.0)
  
    points = sample_fibonacci_sphere_(num_samples)
    for i in range(10): # for the given e1/e2 values 10 iterations are enough
        d = compute_superellipsoid_distance(points[:,0],points[:,1],points[:,2],e1,e2)
        m = np.sum(points**2, axis=1)**0.5 # dist from center
        points = points * ((m-d)/m).reshape(-1,1)
    points *= np.array([a,b,c])
    return points

a,b,c,e1,e2 = 1, 1, 1, 0.5, 0.7
points = sample_superellipsoid_fibonacci_(a,b,c,e1,e2)

fig = plt.figure(figsize=(6,6))
ax = fig.add_subplot(111, projection='3d')
ax.set_title("$a={}$, $b={}$, $c={}$, $e_1={}$, $e_2={}$".format(a,b,c,e1,e2))
_=ax.scatter(points[:,0], points[:,1], points[:, 2], c=-points[:,2], cmap="gist_rainbow")

fig = plt.figure(figsize=(9,9))
e_values = [1.0, 0.5, 0.01]
for i,e1 in enumerate(e_values):
    for j,e2 in enumerate(e_values):
        ax = fig.add_subplot(len(e_values), len(e_values), i*len(e_values)+j+1, projection='3d')
        points = sample_superellipsoid_fibonacci_(1,1,1,e1,e2)
        ax.set_title("$e_1={}$, $e_2={}$".format(e1,e2))
        _=ax.scatter(points[:,0], points[:,1], points[:, 2], c=-points[:,2], cmap="gist_rainbow")
fig.tight_layout()

def sample_superellipsoid_debug_(a, b, c, e1, e2, num_samples=100000, max_iter = 50):
    points = sample_fibonacci_sphere_(num_samples)
    d_ = []
    for i in range(max_iter):
        d = compute_superellipsoid_distance(points[:,0],points[:,1],points[:,2],e1,e2)
        m = np.sum(points**2, axis=1)**0.5
        d_.append(m)
        points = points * ((m-d)/m).reshape(-1,1)
    points *= np.array([a,b,c])
    return points, np.array(d_)

e1,e2 = 0.01, 0.01
points, d = sample_superellipsoid_debug_(1,1,1,e1,e2)
groundtruth = d[-1].reshape(1,-1)
dd = np.max(np.abs(d-groundtruth), axis=1)
plt.xlabel("Iterations")
plt.ylabel("Highest Error")
plt.title("$e_1={}$, $e_2={}$".format(e1,e2))
plt.yscale("log")
plt.grid()
_=plt.plot(dd)

# Optimized version of https://stackoverflow.com/a/26127012
def sample_fibonacci_sphere(samples):
    phi = np.pi * (3 - np.sqrt(5))  # golden angle in radians
    i = np.arange(samples)  
    y = 1 - (i / (samples - 1)) * 2   # y goes from 1 to -1
    radius = np.sqrt(1 - y * y)  # radius at y
    theta = phi * i  # golden angle increment
    x = np.cos(theta) * radius
    z = np.sin(theta) * radius
    return np.array([x, y, z]).T

print("Old fibonacci sphere sampling timings:")
%timeit sample_fibonacci_sphere_(1000)

print("\nOptimized fibonacci sphere sampling timings:")
%timeit sample_fibonacci_sphere(1000)

Old fibonacci sphere sampling timings:
2.85 ms ± 14 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Optimized fibonacci sphere sampling timings:
41.9 µs ± 338 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

def sample_superellipsoid_fibonacci(a, b, c, e1, e2, num_samples=2500, max_iter=10, limit_e1_e2=True):
    if limit_e1_e2:
        e1 = np.clip(e1, 0.01, 1.0)
        e2 = np.clip(e2, 0.01, 1.0)

    phi = np.pi * (3 - np.sqrt(5))        # golden angle in radians
    i = np.arange(num_samples)  
    y = 1 - (i / (num_samples - 1)) * 2   # y goes from 1 to -1
    radius = np.sqrt(1 - y * y)           # radius at y
    theta = phi * i                       # golden angle increment
    x = np.cos(theta) * radius
    z = np.sin(theta) * radius
    s = np.array([x, y, z]).T

    for i in range(max_iter):
        f = (np.abs(s[:,0])**(2./e2) + np.abs(s[:,1])**(2./e2) )**(e2/e1) + np.abs(s[:,2])**(2./e1)
        d = np.sign(f-1) * np.abs(1 - f**(e1/2))
        m = np.linalg.norm(s, axis=1)
        s = s * (1-d/m).reshape(-1,1)

    return s * np.array([a,b,c])

print("Old superellipsoid sampling timings:")
%timeit sample_superellipsoid_fibonacci_(1,1,1,0.5,0.4,1000)

print("\nOptimized superellipsoid sampling timings:")
%timeit sample_superellipsoid_fibonacci(1,1,1,0.5,0.4,1000)

Old superellipsoid sampling timings:
4.27 ms ± 65.1 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Optimized superellipsoid sampling timings:
1.14 ms ± 3.92 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

def sample_superellipsoid_fibonacci_v2(a, b, c, e1, e2, num_samples=2500, max_iter=6, limit_e1_e2=True):
    if limit_e1_e2:
        e1 = np.clip(e1, 0.01, 1.0)
        e2 = np.clip(e2, 0.01, 1.0)

    phi = np.pi * (3 - np.sqrt(5))        # golden angle in radians
    i = np.arange(num_samples)  
    y = 1 - (i / (num_samples - 1)) * 2   # y goes from 1 to -1
    radius = np.sqrt(1 - y * y)           # radius at y
    theta = phi * i                       # golden angle increment
    x = np.cos(theta) * radius
    z = np.sin(theta) * radius
    s = np.array([x, y, z]).T

    for i in range(max_iter):
        f = (np.abs(s[:,0])**(2./e2) + np.abs(s[:,1])**(2./e2) )**(e2/e1) + np.abs(s[:,2])**(2./e1)
        d = np.sign(f-1) * np.abs(1 - f**(e1/2))
        m = np.linalg.norm(s, axis=1)
        s *= (1-(d/m)*(m-d)).reshape(-1,1)

    return s * np.array([a,b,c])

def sample_superellipsoid_v2_debug_(a, b, c, e1, e2, num_samples=100000, max_iter = 50):
    points = sample_fibonacci_sphere_(num_samples)
    d_ = []
    for i in range(max_iter):
        d = compute_superellipsoid_distance(points[:,0],points[:,1],points[:,2],e1,e2)
        m = np.sum(points**2, axis=1)**0.5
        d_.append(m)
        points = points * (1-(d/m)*(m-d)).reshape(-1,1)
    points *= np.array([a,b,c])
    return points, np.array(d_)

e1,e2 = 0.01, 0.01
points, d = sample_superellipsoid_v2_debug_(1,1,1,e1,e2)
groundtruth = d[-1].reshape(1,-1)
dd = np.max(np.abs(d-groundtruth), axis=1)
plt.xlabel("Iterations")
plt.ylabel("Highest Error")
plt.title("$e_1={}$, $e_2={}$".format(e1,e2))
plt.yscale("log")
plt.grid()
_=plt.plot(dd)

print("\nOptimized superellipsoid sampling v2 timings:")
%timeit sample_superellipsoid_fibonacci_v2(1,1,1,0.5,0.4,1000)

Optimized superellipsoid sampling v2 timings:
746 µs ± 12.5 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

from matplotlib.animation import FuncAnimation
from IPython.display import HTML
from scipy import interpolate

def animate(f):
    fig = plt.figure(figsize=(6,6))
    ax = fig.add_subplot(111, projection='3d')

    def update(frame):
        ax.clear()
        f(fig, ax, frame)
        ax.view_init(frame*15+20, frame*60)
        return

    ani = FuncAnimation(fig, update, frames=np.linspace(0, 1.0, 180), init_func=None, blit=False, interval=66)
    plt.close()
    return ani.to_html5_video()

def f(fig, ax, frame):
    n = int(frame*720)
    points = sample_fibonacci_sphere(n)
    ax.set_title("$n={}$".format(n))
    ax.scatter(points[:,0], points[:,1], points[:, 2], c=-points[:,2], cmap="gist_rainbow")

HTML(animate(f))

def animate(f):
    fig = plt.figure(figsize=(6,6))
    ax = fig.add_subplot(111, projection='3d')

    def update(frame):
        ax.clear()
        f(fig, ax, frame)
        ax.view_init(frame*15+20, frame*60)
        return

    ani = FuncAnimation(fig, update, frames=np.linspace(0, 1.0, 180), init_func=None, blit=False, interval=66)
    plt.close()
    return ani.to_html5_video()

def f(fig, ax, frame):
    e1_    = [1.0, 0.5, 0.01, 0.5, 0.01, 1.0]
    e2_    = [1.0, 0.01, 1.0, 0.5, 0.01, 0.5]
    fe1 = interpolate.interp1d(np.linspace(0, 1.0, len(e1_)), e1_)
    fe2 = interpolate.interp1d(np.linspace(0, 1.0, len(e2_)), e2_)
    e1 = fe1(frame)
    e2 = fe2(frame)
    points = sample_superellipsoid_fibonacci_v2(a, b, c, e1, e2, 5000)
    ax.set_title("$e_1={}$, $e_2={}$".format(np.round(e1,2),np.round(e2,2)))
    ax.scatter(points[:,0], points[:,1], points[:, 2], c=-points[:,2], cmap="gist_rainbow")

HTML(animate(f))

Uniform-like Sampling of Superellipsoid Surface

Table of Contents

Parametric Representation¶

The Problem of Sampling with Parametric Representation¶

Solution: Projecting Fibonacci Sphere onto Superellipsoid¶

Fibonacci Sphere¶

Solina's Method¶

Project Points onto Superellipsoid¶

Finding a Good Value for Maximum Number of Iterations¶

Optimized Code: Fibonacci Sphere Sampling¶

Benchmark¶

Optimized Code: Superellipsoid Sampling¶

Benchmark¶

UPDATE v2: Faster and Better!¶

Animated Examples¶

Fibonacci Sphere¶

Superellipsoid¶

Final Thoughts¶

Salih Marangoz